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Abstract 

In common with other iron-based high temperature superconductors, FeSe exhibits a transition to 
a “nematic” phase below 90Kelvin in which the crystal rotation symmetry is spontaneously broken. 
However, the absence of strong low-frequency magnetic fluctuations near or above the transition 
has been interpreted as implying the primacy of orbital ordering. In contrast, we establish that 
quantum fluctuations of spin-1 local moments with strongly frustrated exchange interactions can 
lead to a nematic quantum paramagnetic phase consistent with the observations in FeSe. We show 
that this phase is a fundamental expression of the existence of a Berry’s phase associated with the 
topological defects of a Neel antiferromagnet, in a manner analogous to that which gives rise to 
valence bond crystal order for spin 1/2 systems. We present an exactly solvable model realizing 
the nematic quantum paramagnetic phase, discuss its relation with the spin-1 J\ — J 2 model, and 
construct a held theory of the Landau-forbidden transition between the Neel state and this nematic 
quantum paramagnet. 
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Insulating quantum paramagnets (PM) are magnetic systems with only short-range an- 
tiferromagnetic (AF) correlations even at zero temperature. Examples of quantum para¬ 
magnetic states include valence bond solids (VBS), symmetry protected topological states 
(SPTs), and spin liquids. Interest in quantum paramagnets was greatly intensified following 
the proposal 1 (since shown to be incorrect) that the parent insulator of the cuprate high 
temperature superconductors (HTSC) might be a spin liquid. Because of the proposed rel¬ 
evance to the cuprates, special attention has been paid to square lattice S = 1/2 systems, 
and in particular to the J\ — J 2 Heisenberg model, 

H = J^S.r Sj + J2Y, S i ' S k- (!) 

(h> «jk)) 

where (ij) and (( jk )) denote nearest-neighbor (NN) and second neighbor pairs of sites. 
Numerical studies have shown 2-4 that the ground state of equation (1) has Neel AF order 
for 0 < J2/J1 < 0.4 and stripe AF order [see Fig. 2(a,b)] for 0.6 < Both these phases 

have gapless S = 1 (spin wave) excitations. For 0.4 < J2/J1 < 0.6 the ground state appears 
to have no magnetic order. However, there remain important unresolved issues concerning 
the precise character of this phase or phases; at least for 0.5 < J2/J1 < 0.6 there is fairly 
compelling evidence of translation symmetry breaking VBS order with an energy gap for 
spin-1 excitations. 

The curious fact that upon restoring the spin rotation symmetry the breaking of spatial 
translation symmetry follows rests on a uniquely quantum mechanical effect associated with 
the Berry’s phase of the monopole events. 5,6 A monopole (anti-monopole) is a space-time 
event [see Fig. 1(a)] across which the “skyrmion number” (see the caption of Fig. 1) changes 
by +1(—1). The proliferation of monopoles randomizes the Neel order parameter hence 
causes the system to become a PM. In Ref. 5 , it was argued that monopole events contribute 
to the path integral with a phase factor that depends on the monopole’s spatial location [see 
Fig. l(b,c)]. 


Spin models based on local moments, such as the J\ — J 2 model, have found renewed 
applications in the young field of iron-based superconductivity. The iron-based HTSCs have 
layers of Fe 2+ ions which form a square lattice at high temperatures. Many experimental and 
theoretical studies have concluded that, with the possible exception of heavily phosphorous 
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FIG. 1: The monopole phase factors. A monopole is a singular configuration of the Neel 
order parameter whose direction is depicted by an unit vector h at each point in space and time. 
The Neel order parameter configurations before and after a monopole event differ topologically 

- the skyrmion number changes by one. To understand the skyrmion number imagine assigning 
an unit vector to each point of a two dimensional space subject to periodic boundary conditions. 
Such an assignment is a “map” from a torus to the unit sphere S 2 . These maps can be grouped 
into topological classes, where only maps within the same class can be smoothly deformed into one 
another. An integer, namely the number of times, and sense, the image of the torus wrap around the 
unit sphere, characterizes each class. The skyrmion number is given by q s = -A f dxdy h-d x hxd y n. 
A close inspection of the arrow patterns in panel (a) reveal that before the monopole event (marked 
by the red dot), the skyrmion number is —1, while that after the monopole event it is 0. An 
analogous figure can be drawn for an anti-monopole across which the skyrmion number jumps by 

— 1. Haldane showed that associated with each charge q m monopole event there is a phase factor, 

which enters the path integral over all possible h configurations in space and (imaginary) 
time 5 . This phase factor depends on the spatial location of the monopole core, which it is natural 
to associate with the center of a lattice plaquette, 1Z. There is some arbitrariness in the choice of 
rm, but a consistent pattern for spin-1/2 on a square lattice is shown in panel (b) and for spin-1 
in panel (c). 

doped members of the 122 family, the electronic correlations in the iron-based materials, in 
particular in the iron chalcogenides, are strong ,-9 . Moreover, relatively large local magnetic 
moments in many of these materials have been inferred 10 . For FeSe the itinerant electrons 
only give rise to tiny Fermi pockets 11-14 . Thus, it is plausible that the magnetism of FeSe 
can be addressed using a spin model as a starting point. 

For many iron-based materials, depending on the doping level, there is “stripe” AF order 
at low temperatures [see Fig. 2(a,b)]. Due to the breaking of crystal rotation symmetry the 
stripe order is accompanied by a tetragonal to orthorhombic lattice distortion. In some cases, 
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FIG. 2: The stripe AF and the interactions in H^: Panels (a) and (b) show two degenerate 
versions of a stripe AF ground state corresponding to the two possible directions of the ordering 
vector, where the arrows schematically represent the ordered magnetic moments on the Fe sites. 
Panel (c) illustrates the three-site interaction corresponding to each individual projection operator 
from equation (2). 

e.g. in electron-doped Ba-122 materials, the lattice distortion can exist without stripe order 
[see Fig. 3(a)]. In Ref. 15 it was shown that the distortion is driven by electronic nematicity 
rather than a lattice (phonon) instability. It has been argued that such nematicity reflects 
an underlying stripe ordering tendency, and the reason it can exist without the magnetic 
order is because thermal fluctuations of the continuously varying spin orientation are more 
severe than of the discrete nematic director 16 ~ 18 . 


A potential problem with this perspective is the thermal evolution observed in FeSe 
crystals, in which nematicity onsets at T nem ~90Kelvin, but there is no magnetic ordering 
down to the lowest measured temperatures suggesting the possibility of a zero temperature 
nematic quantum PM phase [Fig. 3(b)]. This fact coupled with the absence of any ob¬ 
served enhancement of the low frequency magnetic fluctuations at T n em 19 ’ 20 has led many to 
concluded that the nematicity in FeSe is driven by orbital ordering. 

Certainly, it is clear that within a local moment picture the lack of magnetic order 
implies that the spin-spin interactions must be highly frustrated; this is true regardless 
of what causes the nematic ordering. Thus in the following we ask, “Can frustrated spin 
interactions alone drive a nematic quantum PM state?” 

We consider models in which each Fe 2+ possesses a localized spin 1 (e.g. the S = 1 version 
of the model in equation (1)). Given the strength of Hund’s coupling and the crystal-field 
splittings, S = 1 is a reasonable possibility for the spin of Fe 2+ ions. However it is important 
to keep in mind that FeSe is not an insulator, and such localized models neglect the effects 
of itinerant carriers. We shall return to this point later. 
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FIG. 3: Schematic phase diagrams: (a) The thermal phase diagram of many iron-based 
superconducting materials. The blue and red dots represent two distinct phase transitions: as a 
function of decreasing T, the discrete crystal rotation symmetry is spontaneously broken at the 
blue point and continuous spin rotation symmetry at the red. (b) A likely thermal phase diagram 
of FeSe. (c) The zero temperature phase diagram of the spin-1 J 1 -J 2 model in equation (1) as a 
function of J 2 /J 1 ■ 


First we demonstrate the existence of a nematic quantum PM phase in an exactly solvable 
Hamiltonian. Consider a square lattice of S' = 1 spins interacting via the short-range, spin 
rotationally invariant Hamiltonian 

H k = KJ2 Ps(Si + Sj + S k ), (2) 

(ijk) 

where K > 0 and Y2{ijk) sums over all elementary triangles of sites [see Fig. 2(c)], and where 
ji and ik are NN bonds and jk is a next-NN bond. Here P 3 (S) = (1/720)S' 2 (S' 2 — 2)(S' 2 — 6) 
is the projection operator onto S = 3. We note that because it involves spin-1 operators, 
equation (2) possesses a global spin SO(3) [rather than SU(2)] symmetry. In addition it 
possesses all crystalline symmetries of the square lattice (i.e. translation and point group 
symmetries). Moreover, there are two degenerate ground states which can be constructed 
exactly as follows: Any closed loop Cqon the lattice can be thought of as a spin-1 chain. 
The famous AKLT (Affleck-Kennedy-Lieb-Tasaki) state 21 of such chains can be written as 
the matrix product state |= El, 1= -i ••• El ; „=-i Tr [M m h)-- A ( m in)} Ki, -,m in ), 
where rn lk = —1,0,1 is the S z quantum number of the spin on site i k , and the matrices 
A(m ) are A(±l) = ( ia y =F l)/2y / 2, A(0) = ct 2 /2, with (Jy : z being the Pauli matrices. We 
identify the two adjacent sites ik and ik+i in an AKLT loop as an “AKLT-entangled pair”, 
and graphically represent it by a blue bond connecting i k and i k +1 in Fig. 4. The two ground- 
states of H k are constructed as the direct product of AKLT loop states on all the loops 
made by connecting nearest-neighbor sites in the x direction, |X) = n^Li |C(i,y),(2 , y ),...,(N, y )), 
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FIG. 4: The nematic paramagnetic ground states of Hk from equation (2). The solid 
blue lines connect AKLT-entangled pairs discussed in the main text. 

or in the y direction, |Y) = Y\^=i \ C (x,i ),(s, 2 ),...,(x,jv)), where (x,y) with x, y = 1,..., N labels 
the sites of a N x N square lattice. The graphical representations of |X) and |Y) are shown 
in Fig. 4. Because the maximum total spin of an AKLT-entangled pair is S — l 21 , it follows 
that the maximum spin component of a triplet of spins (ijk) containing at least one AKLT- 
entangled pair is S = 2. From this it follows trivially that H K annihilates |X) and |Y); 
moreover, since Hk is positive semidehnite, this proves that they are ground states. 


It is only slightly more difficult to see that any other state constructed as a direct product 
of AKLT-loop states has positive energy: since each site can be AKLT entangled with at 
most two other sites, and since there are four elementary triangles per site, in any zero 
energy state exactly two of the sites must be AKLT-entangled, and each entangled pair 
must form an edge of four distinct elementary triangles. The two states |X) and |Y) are 
the only ones that satisfy this constraint. For instance, if the blue bonds in Fig. 4 make a 
turn, an elementary triangle with its vertex at the corner of the turn will not contain any 
AKLT-entangled pair. While we do not have a general proof that no more complex ground 
states exist, explicit finite cluster diagonalization (up to system size 4x4) suggest that there 
are none. 

From the known properties of AKLT states, 21,22 it follows immediately that |X) and |Y) 
are gapped PM states with exponentially falling spin-spin correlations which break crystal 
rotation symmetry: Gx{Rij ) = (X |Si ■ >Sy|X) = 9 cos (Rij ■ QN)8(Rij ■ y)e~^ ij ' x ^R° with 
Qn = 7r (x + y) and £o = 1 / lri(3) in units in which the lattice constant is 1. The gap implies 
that the nematic PM phase is perturbatively stable. However, the asymptotic form of G is 
non-generic - this reflects the fact that Hk lies on a “disorder line 23 ” where, although there 
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is no associated thermodynamic non-analyticity, the oscillatory character of the short-range 
order changes. 

Turning to less “reverse engineered” models, we discuss the very interesting numerical 
(density matrix renormalization group) study by Jiang et al. 2A of the spin-1 version of the 
J\ — J2 model in equation (1). They found that as a function of J2/J1 there is an inter¬ 
mediate paramagnetic phase (characterized by gapped 5=1 excitations) between the Neel 
and stripe ordered phases [see Fig. 3(c)]. Moreover the Neel and stripe order parameters 
vanish continuously as J 2 /Ji approaches the respective critical values from the magnetically 
ordered side (see Fig. 4 of Ref. 24 ). In the supplementary information we perform exact diag- 
onalization on small lattices and demonstrate the adiabatic connectivity of the ground-state 
phases of Hk (equation (2)) and the J\ — J 2 model when J2/J1 falls within the PM regime. 
In addition we compare the excitations in different spin and momentum sectors. 

To obtain an analytic understanding of the nematic quantum PM phase, we consider 
a field theory description valid in the neighborhood of a continuous or weakly first order 
quantum phase transition from the Neel phase to a nematic PM phase. Because the unbroken 
symmetries of the Neel and nematic PM phases do not have a subgroup relationship, classical 
Landau theory would imply that such a continuous phase transition is forbidden. However 
as pointed out in Ref. 25 , when the topological defects of one order carry the quantum number 
of the other order a continuous transition becomes possible. 

According to Ref. 5 , in the path integral describing the quantum fluctuations of the 5 = 1 
Neel antiferromagnet in 2D, the weight associated with each field configuration is determined 
by the usual non-linear sigma model (analogous to the first term in equation (3)), but there 
is also an additional Berry’s phase factor. For a “charge” q m monopole (which causes the 
skyrmion number to jump by q m ) centered on plaquette 1Z (which designates a point on the 
dual lattice) this phase factor (up to a global “gauge” ambiguity) is 77 ^ = e iqrn ^ N ' 7l \ [see 
Fig. 1(c)]. Rotation by 90° about a lattice site transforms 77 ^ —>■ — 77 ^ but keeps q m invariant. 
So long as this symmetry is preserved, the Feynman amplitude of any field configuration with 
an odd q m cancels that associated with the configuration rotated by 90°. As a result odd q m 
monopoles cannot proliferate, although even q m monopoles can. Consequently, states can 
be classified by a conserved skyrmion number parity (—l) 9s . It turns out that the two-fold 
degeneracy associated with opposite skyrmion parity accounts for the two nematic ground 
states depicted in Fig.4. If we introduce a perturbation that breaks the crystal 90° rotation 
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symmetry, the absolute value of the Feynman amplitudes associated with monopoles sitting 
at the +1 and —1 locations in Fig. 1(c) do not need to be the same anymore, hence their 
Feynman amplitudes no longer cancel which means the odd q m monopoles can proliferate 
rendering the PM state non-degenerate. This observation justifies identifying the two-fold 
degeneracy in the symmetric system with spontaneous breaking of 90° rotation symmetry. 

In the following discussions we will use Euclidean space-time and denote the imaginary 
time as r. Consider a 4-component real vector field of unit norm, Cl(x,y,r) = (f2,fhi) with 
|fi | 2 = 1, whose first three components, fl(x,y,r) are the Neel order parameter and 1^4 is 
the Ising-like nematic order parameter. For example, for a spin model on a square lattice, 
we can take H 4 (xi,yi,r) oc (Si ■ S l+ % — Si ■ S i+ y). The non-linear sigma model action we 
shall consider is 


5 = / d 2 xdr 


L2 9 


Id^l 2 + 


+ i 


0 

2ir 2 


d 2 xdr e abcd £} a d x ttbd y Qcd T Qd, 


( 3 ) 


where 0 = ir. This field theory has been introduced in Ref. 25 and Ref. 2 ' which discuss the 
effects of topological terms on the spectrum and phases of non-linear sigma models. For 
our purposes a derivation of equation (3) is given in the supplementary information. Note 
that this particular topological term is possible only in spatial dimension d = 2 with a 
N — 4 component order parameter. Here V is an anisotropy term that favors the first three 
components of fi, (for example V = Afh^] 2 ), hence equation (3) has 0(3) x Z 2 symmetry. 
Consider a monopole configuration centered on the origin, ( x,y,r ) = (0,0,0). Because of 
the anisotropy term, on a two-sphere ( S 2 ) in space-time far away from the center of the 
monopole is largely constrained to lie in the space spanned by the first three components of 
fi, i.e. ~ (h, 0). Topologically the configurations of h on the two-sphere are classified by 
the skyrmion number which, in this case, is equal to the monopole charge q m . However such 
far held configuration can be compatible with two different H 4 orientations in the monopole 
core. Consider the following configurations which reduce to the same q m = 1 configuration 
in h far from monopole the center 


n {±) (x,y,p) = (sin0 ± (p) n(x,y,p),cos0±(p)), (4) 

where p = x 2 + y 2 + r 2 /R, 0+(p) = | tanh(p), and 0_(p) = 7r(l — |tanh(p)). 
Here R is the size of the “monopole core”. both describe q m — 1 monopoles 

but with H = (0, ±1) in the monopole core. It is straightforward to show that 






expjz^fj / cPxdr e abcd Vt a d x Vt b d y Vt c d T Vt d } = e ±J7r / 2 = ±i for these two types of monopole. 
A similar discussion holds for the q m = — 1 monopole. It is easy to generalize the above 
argument to q m = 2 monopoles and obtain the corresponding phase factors e ±l7r = — 1 . 
A consequence of the Z 2 (nematic) symmetry is that the Feynman amplitudes of the two 
monopoles described above have the same absolute values. 

These monopoles represent events (as a function of imaginary time) at which the skyrmion 
number changes. We can thus think of the quantum disordered phase as an interacting fluid 
of skyrmions and anti-skyrmions. Because of the destructive interference between the two 
types of q m = 1 monopoles discussed above, events in which the skyrmion number changes 
by one are forbidden - the net skyrmion number is thus conserved modulo 2 and the Hilbert 
space breaks into an even and an odd sector. Since the stiffness constant 1 /g renormalizes to 
zero in the quantum disordered state and since there is a non-zero (possibly small) density 
of skyrmions and anti-skyrmions, the ground-state energy in these two sectors should be 
equal in the thermodynamic limit. This is the ground-state degeneracy associated with the 
breaking of C 4 symmetry to C 2 . In the supplementary information we discuss the effects of 
an explicit C 4 — > C 2 symmetry breaking field. 

One can also arrive at the PM phase in Fig. 3(c) from the stripe ordered side by prolifer¬ 
ating the monopole of the stripe order parameter. Generalizing the calculation of Ref . 5 we 
find the Berry’s phase factor associated with the stripe monopole is trivial. Consequently 
the charge ±1 monopole can proliferate rendering the resulting quantum disordered state 
non-degenerate. Here, because the symmetries of the two phases have subgroup relationship, 
a continuous transition is allowed in the Landau theory. We expect such transition to be in 
the 0(3) universality class. 

Our key theoretical conclusion is that for a spin-1 model on a square lattice, any “inter¬ 
mediate” quantum disordered phase between the Neel and stripe ordered magnetic phases 
is likely to be a nematic PM. This is consistent with the numerical evidence for the J\ — J 2 
model [equation (1)]. Because the nematic quantum PM phase straddles closely between 
the Neel and nematic phases, we expect it to have low-lying spin excitations near both mo¬ 
menta (7r, 0) and (7T, 7 r). The nematic order considered here might be thought of as “vestigial” 
order 28 left behind when quantum fluctuations have restored spin-rotational symmetry. 

Before speculating on the relevance of these results to the interesting case of FeSe, it is 
necessary to comment on the effects of itinerant carriers. ARPES and STM experiments 
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show very small electron and hole pockets for FeSe 13,14,29 and at the same time quantum 
oscillation associated with these tiny pockets are seen 13,30 . Because the latter experiment 
has very stringent requirement on the long life time of the quasiparticles we take it as an 
indication that the coupling between the itinerant carriers and the local moments is weak. 
This phase with decoupled magnetic and itinerant carrier degrees of freedom is analogous to 
the antiferromagnetic metal phase in heavy fermion systems. However due to the itinerant 
carriers the nematic phase can no longer be characterized as having a spin-gap, and close 
enough to criticality, the universality class of the quantum critical point is likely to be 
altered. In addition, FeSe is at best “quasi-2D,” which is to say that it is ultimately three 
dimensional, and this, too, will alter the nature of the phase transitions, but not necessarily 
affect the sequence of ordered phases. 

Since the iron-based superconductors involve multiple 3 d orbitals it has been long sus¬ 
pected that orbital ordering is the cause of nematicity. Althoug we focus on the nematicity 
caused by frustrated magnetism, we do not intend to imply that orbital degrees of freedom 
play no role at all. As is well known regardless of the driving mechanism the nematic or¬ 
der parameter induces orbital splitting and lattice distortion. Therefore cooperative effect 
involving many different degrees of freedom can well be necessary to describe the system 
behavior quantitatively.In addition the so far ignored spin orbit coupling is ultimately nec¬ 
essary for the nematicity in the magnetic degrees of freedom to induce the band splitting 
observed by ARPES and anisotropic magnetic susceptibility observed by NMR..etc. 

With these caveats, the present results suggest that it may be possible to view the nematic 
phase in the FeSe as being driven primarily by frustrated magnetism. The fact that the 
underlying nematic quantum PM is gapped implies that there need not be any enhancement 
of the low energy magnetic fluctuations associated with the proximity of stripe long range 
order of the sort that has been detected by NMR in other Fe-based materials. Recent 
neutron scattering experiments 31 show relatively short-range magnetic correlations but with 
reasonably large intensity (i.e. with substantial magnitudes of the fluctuating moments) at 
the stripe ordering wavevector at low energies. Consequently, the lack of any strong evidence 
of slow magnetic fluctuations in NMR 20 and the short-range of the magnetic correlations seen 
in neutron scattering appear to be entirely consistent with a magnetic origin of nematicity. 
As discussed previously we expect the nematic PM phase to also exhibit relatively low 
energy spin excitations at the Neel ordering wavevector. The nematic quantum PM state 
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consists of a stacking of spontaneously formed spin-1 AF chains. It can be viewed as a 
“weak” symmetry protected topological state. A potentially testable consequence of this 
observation is that gapless spin 1/2 excitations can arise on certain surfaces and domain 
boundaries. In addition we expect nonmagnetic impurities to cut the chain and induce low 
lying spin excitations in the PM gap. This can be detected by electron spin resonance in a 
manner analogous to that done for spin-1 chains 32 . Further experimental studies of nematic 
order and fluctuations in these materials - including elasto-resistance, nuclear quadrupole 
resonance, and Raman scattering studies - as well more complete neutron scattering studies 
of the magnetic fluctuations would clearly be helpful. 

Finally after the submission of our paper, related theoretical studies appeared 33,34 includ¬ 
ing one based on an itinerant electron picture of the nematicity in FeSe 35 . 
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The Spin-1 Neel-Nematic Transition as a double Spin-1/2 Neel-VBS Transition 

It has been proposed that a Landau-forbidden continuous quantum phase transition be¬ 
tween Neel order and valence bond solid can happen in spin-1/2 systems on square lattice. 
This transition can be described by a non-linear sigma model (NLcrM) with a Wess-Zumino- 
Witten(WZW) term and certain anisotropy terms 36,3 '. The action of this model reads 

Si[(p) = 5 0( 3)xC 4 J0] - 2yri ^7 f dud 2 xdre abcdf (p a d x (j) b dy(j) c d T (j) d d u (j) f . (5) 

Here the 5-component “superspin” </ oc (n x ,n y: n z ,v x ,v y ) consists of the Neel order pa¬ 
rameters n = ( n x ,n y ,n z ) and the two columnar VBS order parameters 
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v = (v x ,v y ) where 


(6a) 

(6b) 


v x ~ ( 1 ) (S(x,y) ■ S(x+l,y) &(x,y) ' &(x—l,y))j 
V y ~ ( — ^) V (S( x ,y) ' &(x,y+ 1 ) — S(x,y) ' S( x ,y-l))- 


u is the auxiliary dimension for dehning the WZW term. So( 3 ) x c 4v is the non-topological 
part of the non-linear sigma model action with 0(3) xC^ v symmetry. This action includes the 
stiffness terms such as f d 2 xdr (^-\dy n \ 2 + ^~\bi addition it contains anisotropy 
terms that favor the Neel order parameters over the VBS order parameters. Moreover, 
among the VBS order parameters there are terms that favor the columnar VBS over the 
plaquette VBS order. Examples of such anisotropy terms include A - (0 2 + 0 2 ) (where A > 0) 
and U ■ (</>! — </>|) 2 (where U < 0). In particular, A > 0 insures that the low energy physics 
is described by the fluctuations of the Neel order parameters. 

Spin-1 can be viewed as two spin-l/2s coupled by strong ferromagnetic(FM) interaction. 
We thus consider two copies of the action equation (5) labeled by superscripts b) and 


sS {1 \ft 2) ] = Si[4> w ]+Si[4> ( - 2) } + 


d 2 xdr (J n n (1) • n^ + J v v ( b . v (2 ^) . 


(7) 


We assume “FM” coupling between the Neel order parameters ( J n < 0) and “AFM” coupling 
between the VBS order parameters (> 0) so that the low energy configurations have 
= n ^ = n and = v. The action in terms of <j> = (n, v) will be similar to 

equation (5) but with a doubled WZW term. 

Note however v = ^(-th 1 ) — v^) cannot be directly measured in spin-1 systems, because 
all physical observables must be symmetric with respect to exchange of the two spin-1/2 
moments. Define two physical order parameters, 


( 1 ) ( 2 ) ( 1 ) ( 2 ) 


v[ = 


Vx Vx 


VyVy 


V v x + 


v' 2 = 


( 1 ) ( 2 ) ( 2 ) ( 1 ) 
Vx J Vy + Vx Vy 


\/ V x+' 


( 8 ) 


v[ carries lattice momentum (0,0), belongs to the B 1 representation of C iv (changes sign 
under 4-fold rotation, but has no sign change under principal axis reflection), and cor¬ 
responds to the nematic order parameter h2 4 defined in the main text (for example, the 
parent Hamiltonian ground state |X) in the main text corresponds to ^ 0 and 

Uyb — — Vy 2 ' = 0 and thus v[ > 0). v' 2 has lattice momentum (7r,7r), is the B 2 represen¬ 
tation of C± v (changes sign under both 4-fold rotation and principal axis reflection around 
a lattice site), and corresponds to certain superpositions of plaquette valence bond solid 
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order (for example, v 2 < 0 may correspond to either v = v y 2 ^ > 0 with plaquette singlets 
centered at (2x + 1/2,2 y + 1/2), or v y 1] = v^ < 0 with plaquette singlets centered at 
(2x — 1/2, 2y — 1/2) for integer x,y). If we parametrize v by (■ v x ,v y ) = (ucos#,usin$), then 
v' = (v[,v 2 ) = (u cos 20, w sin 20). Note that 


v[d y v 2 ~ = 2v 2 d y d = 2 (v x d tl v y - Vyd^v^, 

dyv[ d u v 2 - d y v 2 d u v[ = 2v (d tJ v d u 6 - d tJ 6 d u v ) = 2 (d tJ v x d u v y - d tl v y d v v x ). 


Therefore the WZW term in action S\ in terms of 0' oc (n,v') = (n x ,n y ,n z ,v[,v 2 ) has a 
halved coefficient compared to that in terms of (n,v). The action then becomes 

3 


S'i[0 (1) ,0 (2) ] ~ 5 , i[0 / ] = 5'o(3)xx 2 xz 2 [0 / ]-27ri 


87T 2 


dnd 2 xdr e abcdf (j)' a d x (j)' b d y (j)' c d T ^ d d u (j)' f . (9) 


The action So( 3 )xz 2 xz 2 W\ is derived from 50(3) x c 4 „ [0] while taking equation (8) into account. 
It has 0(3) x Z 2 x Z 2 anisotropy induced by the additional anisotropy terms such as A 4 • 0 4 
and A 5 • 0' 5 2 where A 4! 5 > 0. 

Assume A 5 A 4 > 0. Consider the configuration 


4>{u, x,y,r) = [n(x,y,T)sm(u),cos(u)J , 

where 11 is the 4-component real vector field in Eq. (6) of the main text. When u = 0 this is 
a uniform space-time configuration, and when u — ix j2 this will be a low energy space-time 
configuration with 0g = 0. Integrate over u from 0 to 7 t/2 , the WZW term becomes 

/■tt/2 


— 27T i 
= — 27T i 


d 2 xdr e abcd fl a d x flbd v fl c d T fld 


87T 2 

3 

87T 2 


dw i 


• sin 5 u — sin 3 u cos 2 u) 


d 2 xd t e abcd tt a d x tt b d y tt c d T tt d ■ (—2/3) 


= i 


7T 

2^ 


d 2 xdr e abcd Q a d x Q b d y Q c d r Q d 


In the above the (— sin 5 u) and (— sin 3 u cos 2 u) terms are respectively from abcdf = abcd5 
and abcdf = 5bcda terms in the WZW model (terms with index “5” at other positions 
vanish). This result is exactly the ©-term in equation (6) of the main text with 0 = n. 
Thus the final effective action is given by equation (6) of the main text. (Note that due to 
the anisotropy term V (fl 2 ) the stiffness constant for the first three components of H will be 
different from that of fl 4 at low energies and long wavelengths.) 
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In the above discussion we have made the assumptions that J n < 0 and J v > 0 in 
equation (7), and A 5 A 4 > 0 in equation (9). Here we briefly comment on several other 
possibilities. 

(i) . If J n < 0 and J v < 0, the low energy configurations in equation (7) will be n (1) = 
n^ 2 ' 1 = n and rh 1 ) = v ^ = v, this theory would describe the phase transition from Neel 
AFM order (n) to columnar VBS order ( v ) with a WZW term similar to that of equation (5) 
but with doubled coefficient. 

(ii) . If J n > 0, the low energy configurations in equation (7) would have n ( - 1] = —n^ = n, 
then n is the director of the uniaxial spin-nematic order (uniaxial spin-nematic states have 
|to • S — 0)). This theory would describe the phase transition between ferro-spin-nematic 
order and columnar VBS order (if J v < 0) or a nematic quantum paramagnet (if J v > 0). 
However the WZW term would be absent, so we expect this phase transition to be of first- 
order. 

(iii) . If J n < 0 and J v > 0 as we assumed in equation (7), but A 4 A 5 > 0 in 
equation (9). This theory would describe the phase transition between Neel order and 
plaquette VBS order. 

The effects of explicit C 4 —>• C 2 symmetry breaking field. 

To explicitly break the C 4 symmetry to C 2 , we can introduce a “Zeeman” field 

S —^ S + h J d 2 xdr ^ 4 ( 0 :, y, r), ( 10 ) 

which breaks the degeneracy between the two types of monopole in Eq.(7) of the main text, 
making the absolute value of the Feynman amplitude associated them different - hence they 
do not cancel. Now tunneling events involving unit changes in the skyrmion number are 
allowed, which causes mixing between the even and odd skyrmion sectors which lifts the 
ground-state degeneracy. 

Here we propose the renormalization group flow diagram [Fig. 5] of equation (10). Along 
the vertical axis at h = 0 there is a continuous (Landau-forbidden), or weakly first order, 
phase transition between the small g Neel state and the large g two-fold degenerate nematic 
PM phase. This is supported by Fig. 4 of Ref. 21 . Along the axis where the absolute value 
of h is large the transition from the anisotropic Neel state (where there is Neel long-range 
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FIG. 5: The renormalization group flow diagram of equation (10). The fixed points 
value of g in this figure are all 0(1). The proposed renormalization group flow diagram of 
equation (10). Along the h = 0 axis the red critical point separates the small g Neel ordered state 
from the large g nematic paramagnetic state. Such transition if continuous will be an example of 
Landau-forbidden transition. The blue critical point separates the anisotropic Neel state (a state 
with Neel long range order but (Si ■ Si +X ) + {Si • Si+y)) and the anisotropic paramagnetic phase. 
Such phase transitions should belong to the 0(3) universality class. 

order but (Si ■ S i+X ) ^ (Si ■ S i+y )) into the PM state [this transition happens at the lower 
boundary of the diagonal strip in Fig. 6 discussed below] is in the usual 0(3) universality 
class. In the large g region of Fig. 5, where the system remains PM, tuning h from negative 
to positive encounters a first order phase transition at h — 0 while maintaining a nonzero 
S = 1 gap. 


Returning to the numerical results of Ref. 24 , note that Jiang et al. did study the effect 
of explicit rotation symmetry breaking on their results by introducing anisotropy in the 
NN exchange constants in the x and y directions, so that J\ y > J\ x . They found for 
0 < h = ( J\ y — J\x)/ J\ y < 1 there is always (for some range of J?/■h y ) an intermediate PM 
phase between the Neel and stripe ordered phases [see Fig. 6 for a schematic illustration]. 
Within this PM phase, they found no evidence of a phase transition, suggesting that it 
is all one phase. Remarkably it is found that not only does this PM phase survive for 
0.525 < J 2 /J 1 < 0.555 in the isotropic limit h —> 0, but it also includes the case J 2 = 0 
and h — 1, where the system consists of a decoupled array of spin-1 chains. As it is 
independently known that the spin-1 AF chain is in the same phase as the spin-1 AKLT 
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FIG. 6 : The phase diagram of the J 1 -J 2 model, Eq. 1 of the main text. A schematic 
reproduction of the phase diagram reported in Ref. 24 . The line of J\ y — J\ x = 0 shows the phase 
diagram for the J 1 -J 2 model, Eq. 1 of the main text, with fourfold rotation symmetry. 

chain, this observation nicely connects the results of the J\ — J 2 model to those obtained for 

h k . 


Finally, the other interesting scenario 2, , namely, the gapless state of Eq. (6) in main text 
at 0 = 7T can be viewed as the critical state between two different SPTs (one at 0 = 0 
and the other 0 = 2n). Presumably such critical state can be obtained by proliferating the 
domain walls in the nematic order parameter in our nematic PM phase. Further works are 
definitely deserved. 

Exact diagonalization results for the the 5=1 J 1 -J 2 model 

Exact diagonalization was performed on 4 x 4 square lattices with periodic boundary 
conditions for the spin-1 Ji-J 2 antiferromagnetic Heisenberg model 

H = S > • Si + J 2 J2 S -i ■ S '*• (U) 

{*?') <(i fc >> 

The spin singlet and spin triplet gaps at momenta q = ( 7 r, 0) and (7r, 7 t) are presented for 
different values of J 2 /Ji in Fig. 7 The global spin-1 gap is given by the minimum of the 
blue and dashed black curves. As a result it exhibits a kink consistent with that reported in 
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Ref. 21 . What is noteworthy are (1) the spin-0 gap plunges in the range of J 2 /J 1 where the 
nematic quantum PM state is expected. This presumably reflects the small splitting (due to 
quantum tunneling) between the two states that would be degenerate in the thermdynamic 
limit. (2) As the (7r, 0) triplet gap vanishes as J 2 /J 1 approaches the PM to stripe phase 
boundary, the (vr, 7r) triplet gap steadily increases. The reverse is true as J 2 /J 1 approaches 
the PM to Neel phase boundary. That within the PM regime, the A,g = i(7r,7r) is small 
compared to As=i(7 T, 0) on the small J 2 /J 1 and large on the large J 2 /J 1 side suggests the 
existence of two closeby quantum phase transitions plays a key role in the physics. 

The ground state fidelity susceptibility 38 is presented in Fig. 8. This quantity displays a 
clear peak providing strong evidence of a quantum phase transition(s) within 0.5 < J 2 /J 1 < 
0.6, even on such a small lattice. This result could be interpreted as indicating a single 
strongly first-order transition (which is inconsistent with the result of Ref. 24 ). Alternatively 
it can be taken as evidence of the existence of two (continuous) transitions (favored by the 
result of Ref. 24 ). 


Exact diagonalization results for a model interpolating between the parent Hamil¬ 
tonian and J 1 -J 2 model 


We consider a model which interpolates between the parent Hamiltonian Hk in Eq. (2) 
of the main text and the J\-J 2 Heisenberg model, 

15 Ji 


H\=\- 


4 A 


■H k + ( 1 -AMA + 52 Si-Si) 


(h> 


««» 


Jx 52 S% ' s 3 + [ J 2 + a (Ji/2 - J 2 )\ Si ■ Sj T A • (higher order terms), 

(ij) «h» 


( 12 ) 


where the “higher order terms” contain those terms involving 4 or 6 spin operators. This 
model becomes the parent Hamiltonian AAA • H K at A = 1, and the J\-J 2 Heisenberg model 
at A = 0. We study the behavior of this model for three different J 2 /J 1 values, 

• J 2 /J 1 — 0: the J 1 -J 2 Heisenberg model in this case should exhibit Neel order. The 
results for the interpolating model are shown in Fig. 9 and Fig. 10. In particular the 
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FIG. 7: The global singlet(As' = o) gap and the triplet(As=i) gaps at momenta q = (7 r, 0) and (n, 7 r) 
for the spin-1 J 1 -J 2 model on 4 x 4 lattice obtained by exact diagonalization. The global triplet gap 
result is the minimum of the blue and black dashed curves. It exhibit a sharp kink consistent with 
the DMRG results of Ref. 21 . Due to the quantum tunneling between the two degenerate nematic 
PM states on finite lattices, one expects an unique singlet ground state and a small gap for the 
singlet excitations. The unique ground state for a given a = J 2 /J 1 is a spin singlet and has lattice 
momentum (0,0). The singlet gap is small when J 2 /J 1 falls in the region where the nematic PM 
state exists (it nearly vanishes around a = 0.6 where A 5 = 0 ( 0 ; = 0.6) = 0.0033Ji). Moreover the 
lowest energy singlet excited state has lattice momentum ( 0 , 0 ), suggesting no translation symmetry 
breaking in the tentative nematic quantum PM state. The triplet gap has a cusp around a = 0.54. 
The lowest energy S = 1 states for a < 0.54 have lattice momentum (it, 7r) consistent with Neel 
order, and the lowest energy S = 1 states for a > 0.54 have lattice momentum ( 7 r, 0) or (0, ir) 
consistent with stripe antiferromagnetic order. 

ground state fidelity susceptibility shown in Fig. 10 has a prominent peak at around 
A = 0.9, suggesting a phase transition from nematic paramagnet at A = 1 to Neel 
order at A = 0. 


• J 2 I — 0.54: according to the DMRG result in Ref. 24 , the J 1 -J 2 Heisenberg model in 
this case will be in PM phase. The results for the interpolating model are shown in 
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FIG. 8: Ground state fidelity susceptibility for the spin-1 J 1 -J 2 model on 4x4 lattice. The ground 
state fidelity Fo(a, a + 5a) = |(^o(a: + 5a)|^o(a0)| is the overlap between ground states (V’o) at 
parameter a and a + 5a. There is a sharp peak at around a = 0.56. 



FIG. 9: Singlet(As = o) and triplet(Ag = i) gaps for the spin-1 interpolating model equation (12) 
with J 2 = 0 on 4 x 4 lattice. 

Fig. 11 and Fig. 12. The ground state fidelity susceptibility shown in Fig. 12 has no 
peak, suggesting that the nonmagnetic phase at A = 0 (J 2 /J 1 = 0.54) is also a nematic 
paramagnet as the ground states of A = 1 parent Hamiltonian. 
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FIG. 10: Ground state fidelity susceptibility for the spin-1 interpolating model equation (12) with 
J 2 = 0 on 4 x 4 lattice. There is a peak at around A = 0.9, suggesting that this marks the transition 
point of a nematic paramagnetic phase for A > 0.9 to a Neel ordered state for A < 0.9. 
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FIG. 11: Singlet(As=o) and triplet(As = i) gaps for the spin-1 interpolating model equation (12) 
with J 2 /J 1 = 0.54 on 4 x 4 lattice. The singlet gap remains small, suggesting that the nematic 
paramagnet phase persists to J 1 -J 2 model limit (A = 0). 

• J 2 I Ji = 1: for this parameter choice the J 1 -J 2 Heisenberg model should be deep inside 
the stripe magnetic ordered phase. The results for the interpolating model are shown 
in Fig. 13 and Fig. 14. The ground state fidelity susceptibility shown in Fig. 14 has 
a prominent peak at around A = 0.8, suggesting a phase transition from nematic 
paramagnet at A = 1 to stripe order at A = 0. 
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FIG. 12: Ground state fidelity susceptibility for the spin-1 interpolating model equation (12) with 
J 2 /J 1 = 0.54 on 4 x 4 lattice. There is no peak, suggesting no phase transition from nematic 
paramagnet at A = 1 to the nonmagnetic phase at A = 0. 



FIG. 13: Singlet(As=o) and triplet(A,s = i) gaps for the spin-1 interpolating model equation (12) 
with J 2 /J 1 = 1 on 4 x 4 lattice. The singlet gap remains small for A > 0.8, suggesting that the 
nematic paramagnet phase persists in this region. 
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FIG. 14: Ground state fidelity susceptibility for the spin-1 interpolating model equation (12) with 
J 2 /J 1 = 1 on 4 x 4 lattice. There is a peak at around A = 0.8, that this marks the transition point 
of a nematic paramagnetic phase for A > 0.8 to a stripe ordered state for A < 0.8. 
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